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1.  Introduction 


In  physical  geodesy,  the  determination  of  various  quantities,  like  the  geoid  undulation,  the 
deflection  of  the  vertical,  the  terrain  correction,  and  similar  quantities  attributable  to  mass  attraction, 
involve  the  calculation  of  convolution  integrals.  In  these  applications,  one  may  call  the  two 
functions  being  convolved  the  signal  function  {or  data)  and  the  kernel  function,  where  the  latter 
represents  the  model,  or  system,  that  transforms  the  data  in  some  physically  meaningful  way.  For 
example,  Stokes’  formula  is  a  convolution  of  gravity  anomalies  (the  data)  with  Stokes’  function 
(the  kernel)  that  produces  the  geoid  undulation.  Although  these  convolutions  are  of  functions  on 
the  sphere,  only  planar  approximations  will  be  considered  here. 

It  is  well  known  that,  because  the  data  that  enter  the  convolution  are  finite  and  discrete  and 
because  certain  fast  computational  techniques  require  it,  the  convolution  theoretically  defined  on  the 
entire  plane  must  be  replaced  with  a  convolution  for  periodic  functions.  This  use  of  the  so-called 
circular  convolution,  instead  of  the  linear  convolution,  therefore,  naturally  introduces 
corresponding  errors  in  the  model  for  the  desired  quantity. 

There  are  two  basic  approaches  to  treating  this  error  that  have  been  used  in  the  past.  The  first 
simply  ignores  the  error.  This  was  the  case,  for  example,  in  Schwarz  et  al.  (1990);  and  it  has 
justification,  as  will  be  seen.  The  second  approach  is  based  on  extending  the  finite,  discrete  data 
arrays  being  convolved  with  zeros.  This  artifice  was  discussed  by  Oppenheim  and  Schafer 
(1975).  Many  numerical  tests  have  been  done  to  show  that  this  so-called  zero  padding  improves 
the  computation  of  Stokes'  integral  in  planar  approximation  using  the  discrete  Fourier  transform 
(or,  equivalently,  the  fast  Fourier  transform,  FFT);  one  may  cite  Tziavos  (1992),  Sideris  and  Li 
(1993),  and  Haagmans  ct  al.  (1993).  Furthermore,  several  padding  variations  have  been 
investigated  numerically,  including  the  extension  (padding)  of  both  the  data  and  kernel  with  zeros 
(Zhang,  1995),  the  extension  of  the  data  with  zeros  and  the  kernel  with  its  known  values  (Sideris 
and  Li,  1993),  as  well  as  the  tapering  of  the  data  in  a  border  to  zero  at  the  edges  (Vermeer,  1995; 
Tziavos,  1996).  It  is  clear  from  this  recent  geodetic  literature  that  there  remains  some  confusion  as 
to  the  nature  and  appropriateness  of  the  padding  schemes  and  the  relationship  between  the  linear 
and  circular  convolutions.  It  is  noted  that  Haagmans  et  al.  (1993)  make  the  clearest  statement  that 
the  correct  padding  scheme  is  one  where  only  the  data  are  extended  with  zeros  and  the  kernel  is 
naturally  extended.  However,  they  give  only  a  graphical  “proof’. 

The  purpose  of  this  report  is  to  analyze  the  padding  schemes  that  do  not  yield  equality  between 
the  circular  and  linear  convolutions,  and  to  determine  the  quantitative  relationship  between  the 
circular  convolution  error  and  the  edge  effect.  Explicit  formulas  for  the  errors  in  the  padding 
schemes  and  for  the  edge  effect  are  derived;  and  these  errors  are  then  analyzed  using  transfer 
functions  and  a  numerical  example.  The  analysis  proceeds  to  properly  defined  frequency 
responses  of  discrete  derivative  operators,  including  those  that  are  applied  to  convolutions.  The 


mathematical  derivations  are  restricted  to  a  single  dimension  in  order  to  maximize  simplicity.  The 
results  can  easily  be  generalized  by  inference  to  higher  dimensions  and  are  given  explicitly  for  two 
dimensions. 


2.  Convolutions 


Three  types  of  convolution  will  be  needed.  The  first  is  for  continuous  functions  defined  over 
the  real  line.  Let  g(t)  and  h(t)  be  two  functions,  where  t  is  a  real  number.  Their  (linear) 
convolution  is  denoted  by  (g  *  h)(t)  and  is  given  by 


oo 

r 


(g*h)(t)  = 


g(t)  h(t  -  T)  dT  , 


J 

— CX5 


— OO  <  t  <  oo 


(1) 


It  is  assumed  that,  as  1 1 1  ^  ^ ,  the  functions  g(t)  and  h(t)  attenuate  to  zero  in  such  a  way  that  the 
convolution  (1)  exists  for  all  t. 

For  discrete,  infinite  sequences,  g|^  and  h|^ ,  of  presumably  evenly  spaced  samples  of  these 
functions,  the  linear,  discrete  convolution  is  denoted  by  (g  #  h)|^ ,  and  it  is  given  by 

oo 

(g#h)k=  X  gnhk_„  ;  ke  Z  (2) 

n 

where  Z  is  the  set  of  integers;  and,  again,  it  is  assumed  that  the  sum  exists. 

Finally,  consider  the  periodic  sequences,  (gNj^  having  the  same  period,  N. 

The  convolution  as  defined  in  (2)  does  not  exist  for  these  sequences,  and  an  alternative  definition  is 
required.  The  periodic  (or  cyclic),  discrete  convolution  is  denoted  by  ^gN#hN)k’ 
given  by 

N/2-  I 

(iNShN),'  X  (iN)„(iiNL,. ;  -f<k<f-,  ,3) 

n  =  -N/2 

where,  without  loss  of  generality,  one  may  assume  that  the  integer  N  is  even.  Any  consecutive  set 
of  integers  n  in  the  summation  may  be  used  because  of  the  periodicity  of  the  sequences.  Of 
course,  the  convolution  itself  is  also  periodic  with  the  same  period,  N. 

The  discrete  Fourier  transform  (DFT)  (gwjk  defined  by 


-2- 


(4) 


(GN),'DFT(8N)-”|j(iN)„e"T“;  (.0,...,N-1 

Here,  the  more  traditional  index  sequence,  k,  (•  =  0, ...,  N  -  1  ,  is  used;  however,  because  of  the 
periodicity  (also,  of  the  exponentials),  any  other  sequence  of  N  indices  yields  an  equivalent 
transform.  Because  (g^)^  is  assumed  to  be  a  real-valued  sequences,  the  DFT  has  the  following 
conjugate-symmetry  property,  that  can  be  proved  easily  from  (4): 

(gN)i^is/-ea/forallk  o  (g =  (Gn)^_j  foranyC  (5) 

where  superscript-*  signifies  complex  conjugate. 


3.  Approximations 


The  continuous  convolution  (1)  represents  the  reality  (the  model)  of  the  operation  that 
combines  the  data  and  kernel  functions,  g  and  h.  In  practice,  one  has  a  finite  number  of  discrete 
values  of  the  signal,  that  is,  a  finite  sequence  of  N  (equally  spaced)  data,  g|^ .  If  these  data  are 
given  (again,  without  loss  in  generality)  on  the  interval  -N/2<k<N/2-  1  ,  then  the  true 
convolution  (1)  may  be  broken  into  a  number  of  parts: 


T/2 


g(x)  h(t  -  X)  dx  ^ 


g(x)  h(t  -  X)  dx  + 


g(x)  h(t  -  x)  dx  + 


g(x)  h(t  -  X)  dx 


-T/2 


T/2<|Ti<T/2+]t| 


T/2+|t|<|t|<oo 


(St  h)(0  ^cdgc  ^trunc 


(6) 


(Sn  ^  ^cdge  ^discrct  ^tRinc 


where  T  is  the  corresponding  continuous  interval  on  the  real  line  (if  At  is  the  sampling  interval, 
then  T  =  N  At ),  and  where  the.  function,  g-|-(t) ,  is  given  by 


g(t), 

0 ,  otherwise 


(7) 


The  truncated  sequence,  (g^jk  ,  represents  the  collection  of  available  samples  of  the  signal,  i.e., 
the  data: 


-3- 


(8) 


(gwlk 


gk. 

0, 


_N<k<N_ 
2  2 

otherwise 


The  first  part  of  the  final  equation  in  (6),  the  discrete  convolution,  can  be  calculated  from  the  data, 
while  the  second,  third,  and  fourth  parts.  Sedge  >  ^discrcf  ^trunc’  represent  errors,  called, 
respectively,  the  edge  effect,  the  discretization  error,  and  the  truncation  error  (see  Figure  1). 


Figure  1:  Truncation  error  and  edge  effect  for  the  convolution  (gj  *  h)(t) .  Also  shown 
for  later  reference  is  the  cyclic  convolution  error. 

The  edge  effect  and  truncation  error  are  due  to  the  finite  extent  of  the  data,  but  the  edge  effect  is 
considered  separately  since  it  can  be  avoided  to  some  extent  by  limiting  the  domain  of  t  to  an 
interval  smaller  than  the  data  interval  (it  is  zero,  when  t  =  0 ).  Neither  the  discretization  error  nor 
the  remaining  truncation  error  is  within  the  scope  of  this  discussion.  Both  can  be  reduced  to 
acceptable  levels  with  sufficient  sampling  in  density  and  extent.  Instead,  the  error  in  computing  the 
convolution  by  DFT  is  now  investigated.  It  is  assumed  that  the  kernel  function,  h,  is  known  for  all 
t;  therefore,  it  can  be  sampled  for  all  k. 

In  order  to  use  the  DFT  (or,  FFT),  one  must  assume  that  both  the  data  and  the  kernel 
sequences  are  periodic.  Thus,  to  proceed,  consider  the  following  assignments: 


(iwlk^gk^ 

(gN)k  +  l'N  Sk  ’ 

-^<k<N-_i  , 

I'e  Z 

(9) 

nIk  ~  k  ’ 

(^N)k  +  (N  ^k  ’ 

-^<k<N_l 

,  I'e  Z 

(10) 

Then,  the  linear,  discrete  convolution  in  (6)  should  be  approximated  by  a  cyclic,  discrete 
convolution.  There  are  the  following  alternatives,  in  each  for  -  N/2  <  k  <  N/2  -  1  : 


-4- 


(11) 


(gN^h)k“(8N^hN)|^-(ei)i^ 

(gN^h)k  =  (g2N#h2N)k-(e2)k  (12) 

(gN#h)k  =  (g2N^h2N),^-(e3)|j  (13) 

where  g^'fg  (similarly,  hj^)  is  a  periodic  sequence  defined  by  the  extension  of  g,sj  with  zeros  on 
either  side  for  an  interval  equal  to  half  the  length  of  the  original  sequence: 


gk>  -y^k<^-i 
0,  -N<k<-^-l  and  ^<k<N-l 


(g2N)k  +  2fN  “  (g2N)k  ’  -N<k<N-l  ,  PeZ 


(14) 


Also,  the  extended  periodic  sequence  h2N  is  defined  by 

(h2N)k  =  hk’  (h2N)k  +  2eN^(^2N)k  ’  -N<k<N-l,  PeZ  (15) 

The  first  approximation  (11)  appears  in  the  early  geodetic  literature  dealing  with  Fourier 
techniques  (e.g.,  Schwarz  et  al.,  1990).  Equation  (12)  was  considered  by  some  as  a  way  to 
reduce  the  circular  convolution  effects  through  zero-padding  of  both  data  and  kernel  functions. 
The  third  possibility,  (13),  reflects  the  case  where  the  data  are  extended  with  zeros,  but  the  kernel 
sequence,  though  assumed  periodic  (period  2N),  is  extended  naturally  using  its  own  values. 
Elaborating  on  this  case,  it  is  easy  to  prove  that  £3=0.  First  note  that 

(g2N;n  =  lgN)n  ;  -N  <  n  <  N  -  1  (16) 

(h2N)k  .  n  =  (h2N)k-n  =  hk-n  ;  -  ^  <  H  <  ^  -  1  and  -  ^  <  k  <  ^  -  1  (17) 


-5  - 


Therefore  we  have,  for  -  N/2  <  k  <  N/2  -  1  : 


(g2N#h2N)k=  2  (gN)n(h2N)k-n  ^  from  (3)  and  (16) 

n  =  -N 

N/2-  1 

=  E  (gN)nhk-n  ;  from  (8)  and  (17) 
n  =  -N/2 

Hence,  from  (8)  and  the  definition  of  discrete  convolution  (2), 

(gN^h)|^  =  (g2N#h2N)k  ! 

showing  that  £3  =  0,  i.e.,  that  the  cyclic  convolution  of  the  zero-padded  data  and  the  naturally 
extended  kernel  sequences  equals  the  linear  convolution  of  the  original  sequences. 

This  result  can  be  specialized  immediately  to  the  case  where  the  kernel  sequence  is  also  padded 
with  zeros  (like  g2N  in  (14)),  because  the  natural  extension  of  h^  is  h2N  • 

(gN^hN)k  =  (g2N^h2N),j  ;  -y<k<y-l  (19) 

Both  equations  (18)  and  (19)  relate  linear  convolutions  to  corresponding  cyclic  convolutions. 
Equation  (19)  is  the  justification,  originating  in  Oppenheim  and  Schafer  (1975),  for  padding  both 
sequences  with  zeros  over  appended  intervals  on  either  side  equal  in  length  to  N/2.  But  this  is  the 
correct  procedure  only  if  both  sequences  are  available  only  over  the  truncated  interval 
-  N/2  <  k  <  N/2  -  I  .  It  is  clear  from  (6)  that,  since  we  know  the  function  h  (it  need  not  be 
truncated),  the  proper  cyclic  convolution  to  be  used  is  (18),  not{\9). 

It  can  be  shown  with  due  consideration  of  the  sequence  definitions,  that  the  errors  associated 
with  the  alternatives  (1 1)  and  (12)  are  given,  respectively,  by 


1  —  1 

E 

gn 

k  -  n  +  N 

^k-n)  ’ 

-^<k<-2 

n  =  k+N/2+l 

0 

k  =  -l 

-N/2  +  k 

0<k<^ 

1  E 

gn 

(hk-n-N 

-hk-n)' 

n  =  -N/2 


and 


-6- 


/ 


(e2)k=\ 


N/2  -  1 

2  gnhk-n> 

n  =  k+N/2+l 

0 

-N/2  +  k 

“  ^  gphk-n’ 

n  =  -N/2 


-^<k<-2 
k  =  -l 

0<k<^-  1 


(2i) 


The  fact  that  the  errors  in  these  cases  are  zero  at  k  =  -1  instead  of  k  =  0  is  an  artifact  of  N  being 
even.  Note  that  the  summand  in  (6|)|^  (see  Figure  1)  involves  the  values  of  the  kernel  function 
potentially  close  to  the  origin,  while  for  ,  only  the  tail  ends  of  h  enter.  Hence,  if  the  kernel  is 
largest  near  the  origin  and  attenuates  with  distance  from  the  origin  (many  geodetic  kernels,  in  fact, 
behave  like  the  reciprocal  distance),  then  the  error  (e  |  )k  is  generally  much  larger  than  (E2)k  •  Of 
course,  this  comparison  is  almost  irrelevant  if  alternative  (18)  is  used.  However,  each  method 
differs  in  the  computational  load  and  computer  memoiy  requirements,  and  these  may  be  important 
considerations. 


4.  Edge  Effect 


The  edge  effect  error,  for  0  <  t  <  T/2  ,  is  given  by 


T/2+t 


^cdge^O 


g(T)  h(t  -  X)  dx 


T/2 


(22) 


with  a  similar  integral  for  -  T/2  <  t  <  0 .  An  approximation  is  a  discretization  of  this,  again,  for 
0<k<N/2-l  : 


N/2  +  k 

(^cdge)k  ^  Sn^k-n  (23) 

n  =  N/2 

and  a  similar  expression  for  -N/2  <  k  <  -2  .  The  edge  effect  is  zero  for  k  =  -1  .  It  is  already  clear 
from  Figure  1  that  (ecjggji^  and  (e|)k  have  commensurate  error  characteristics. 

Figure  2  shows  that  the  edge  effect  and  truncation  error  remain  unchanged  with  the 
introduction  of  the  extended  sequence  g2N .  The  cyclic  convolution  error  in  this  case  is  zero  (as 
shown  in  ( 1 8))  because  the  extended,  periodic  kernel  function  multiplies  only  zeros  in  the  extended 
data  sequence. 
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Figure  2:  Truncation  error,  edge  effect,  and  cyclic  convolution  error  for  the  extended, 
zero-padded  data  sequence. 


5.  Extension  to  Two  Dimensions 


The  derivations  of  the  previous  section  for  one-dimensional  functions  and  sequences  can  easily 
be  extended  to  two  dimensions,  if  the  coordinates  are  Cartesian.  The  appropriate  cyclic 
convolution  to  use  in  place  of  the  linear  convolution,  from  ( 1 8),  is  formally  given  by 


^  -  (gS„2N,  #  i 


N|  ,  N, 


iij<k 
2  -^2-  2 


where 


(24) 


-(),()  \ 
t?2N„:N,L 


NJ 

0. 


and  -  ^  <  k  2  ^  ■ 


N, 

~T 


N2<k2<-^-l  or^^<k2<N2-l 


I'D.O  ] 

lg:N,.2NJt,.2f,N„k2  +  2t2N2~l^2N, 


,2N;) 


k  i.k. 


■N|<k|<N|-l,  -N2<k2<N2-l,  ('i.lS^Z 


(25) 


That  i.s.  the  zero-padded  signal  array,  ,2n,  ’  the  original  array  plus  a  border  of  zeros,  whose 
width  is  either  (N|l/2  or  (N2)/2,  depending  on  the  coordinate  direction.  This  extended  array  is 
continued  periodically  over  the  entire  plane.  The  periodic  kernel  array  is  defined  analogous  to  ( 1 5): 
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-N,  <k,  <N|  -  1  ,  -N2<k2<N2-  1 


(h2N„2N3)k„kr^'''-''^’ 

(h2N„2NJk,+2(,N„k,  +  2t,Nr(^2N,.2NJk„k,’ 

(26) 

where  the  extension  to  the  larger  2N  |X2N2  grid  is  accomplished  using  the  actual  values  of  the 
kernel  function. 


6.  Alternative  Indexing 


Most  FFT  algorithms  assume  the  DFT  is  defined  with  indices  starting  at  zero,  as  in  (4). 
Because  of  its  periodicity,  the  cyclic  convolution  that  is  identical  to  the  linear  convolution  still  is 
Eq.  (23).  The  difference  is  in  the  padding  of  the  extended  kernel  array  prior  to  convolution.  In 
essence,  the  padding  scheme  is  no  different  than  what  is  specified  by  (26),  but  by  shifting  the 
index  to  start  at  zero,  the  extended  part  of  the  array  is  not  filled  with  the  natural  values  of  h,  but  by 
the  values  of  h  for  negative  indices.  That  is,  we  simply  use  an  interval  of  the  domain  of  h2N ,  2N2 
other  than  the  principal  one  given  in  (26),  that  happens  to  start  at  (0,0). 

One  has  the  following  algorithm  for  padding  the  data  and  the  kernel,  respectively: 


[  Ski.kj  ’ 


0<k,<N,-l,  0<k2<N2-l 

N,<k,<2N,-l,  0<k2<N2-l 
0<k|<N,-l,  N2<k2<2N2-l 

N,<k,<2N,-l,  N2<k2<2N2-I 


(27) 


k  i,k2  ’ 


(^2N,.2nJ 


2N,.2N2jk„k2"\  h 


'k|-2N„k2  > 

k„k,~2N,  ’ 


'k|-2N|,k2-2N2  ’ 


0<k,  <N,  -  1  , 

N,  <k,  <2N,  -  I  , 
0<k|  <N,-1  , 

N|  <k|  <2N,  -  I  , 


0<k2<N2-  1 
0<k2<N2-  1 
N2<k2<2N2- 1 
N2<k2<2N2- 1 


(28) 


7.  Error  Analysis  Equations 

While  it  is  easy  to  justify  mathematically  that  (24)  should  always  be  used  to  compute  the  linear 
convolution  (with  zero  error),  the  tax  on  computer  memory  and  time  may  be  prohibitive  for  very 
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large  arrays  (although,  computer  memory  and  processing  speeds  are  still  increasing  at  near 
exponential  rates).  For  this  reason  it  is  instructive  to  study  the  errors  associated  with  the  use  of  the 
approximation  corresponding  to  (1 1): 


(gN,.N^#h)  -(gN,,N 


-.2 


k  ,,k 


-(ei) 


k  i»k2 


2  Sk,  s  2  » 


(29) 


And,  though  it  yields  no  savings  in  computer  memory,  the  use  of  zero-padded  kernels  does  offer 
some  savings  in  computational  time  (for  dimensions  greater  than  1,  only),  since  the  computation  of 
DFT’s  of  zero-vectors  in  the  padded  array  is  trivial.  Thus,  also  the  second  approximation  may  be 
considered,  corresponding  to  (12); 


(gN,.N2#h' 


k  i.kj 


-0.0 

gN„N 


“N  N. 


k,.k2 


(e 


2ik„k2 


-V<k 


4-. 


N, 


<k,< 


(30) 


In  the  following,  the  errors  (ei),^  and  (£2)^,  kj  ’  effect,  are  characterized  in 

terms  of  their  transfer  functions  -  how  they  affect  the  convolution  in  terms  of  signal  frequency. 

Explicit  expressions  for  these  errors  are  analogous  to  equations  (20),  (21),  and  (23). 
However,  in  two  dimensions,  nine,  instead  of  three,  different  expression  are  obtained,  depending 
on  the  domain  of  (k  |,k2) .  In  addition,  there  are  up  to  three  summations  over  different  ranges  for 
each  subdomain  of  (k  ,,k2) .  To  simplify  the  notation,  the  following  generalized  expression  is 
used: 


•^2  ni^2 


(31) 


where  the  sum,  S ,  represents  a  collection  of  sums,  each  over  two  indices  of  possibly  different 
ranges;  and  where  Uk  ^ kj  -  n2  depends  on  values  of  the  kernel  function,  as  follows: 


(^1  Iki.kj- 

(“')k|-n„k2-n2^ 

(32) 

(^2)k„k2- 

(‘^2)k,  ..n,.k2  - 

(33) 

(Ecdgclk^k,' 

(Ucdgc)|,^_  n„k2  ■  "'■'^2  "2 

(34) 

In  Table  1,  the  ranges  of  (n ,,  n2)  are  listed  for  a  partial  set  of  the  possible  domains  for  (k,,k2)  in 
the  case  of  cyclic  convolution  errors,  Eqs.  (32)  and  (33).  Similar  index  ranges  can  be  derived  for 
the  other  domains  of  (k  |,k2) ,  but  one  may  limit  the  further  discussions  to  those  listed  if  the  kernel 
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function  is  symmetric  or  antisymmetric  (as  it  is  in  geodetic  applications).  Table  2  lists 
corresponding  ranges  of  indices  for  the  edge  effect  (34). 

Table  1;  Ranges  of  indices,  (nj,  n2) ,  and  values  of  a^,  a2  (in  (32)),  for  sums  of  (32)  and  (33) 
entering  in  equation  (31).  Domains  -  N|/2  <  k ,  <  -  2  and  -  N2/2  <  kj  <  -  2  are  omitted. 


It  is  now  assumed  that  the  signal  is  a  stationary,  stochastic  process  on  the  plane  with  zero 
mean.  Using  the  notation  in  (31),  the  variance  of  the  error  E]^  j.^  is  then  given  by 


^  ^  ^k,-n|,k2-n,  ^ 


(35) 


n,,n2  * 


fen  |,n,  fen  n. 


where  E[-]  is  the  statistical  expectation  operator. 
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Table  2:  Ranges  of  indices,  (n|,  02) ,  for  sums  of  (34)  entering  in  equation  (31).  Domains 
-  N ,/2  <  k  I  <  -  2  and  -  N2/2  <  kj  <  -  2  are  omitted. 


E  gn,,n2  8n'  n'  auto-covariance  of  the  discrete  signal,  which  is  also  the  continuous- 

signal  auto-covariance  function  sampled  at  the  discrete  points  (n  |  -  n  ,,n2  -  n2) .  If  (t)g(fi  ,f2)  is  the 
power  spectral  density  of  g,  depending  on  spatial  frequencies,  f|,f2 ,  then 


E[gn.,n2  8„;,n;]=  (D,(f„f2)  e' "^>1  df,  df2 


coordinates  (k  |,k2) . 

With  suitable  reorganization  of  indices,  it  is  easily  shown  that 
Uk„k,(fl.f2)  u;„i/f„f2)  =|Uk„|„(f„f2)p 

=  2  u_  „  cos 

k,-n„k2-n2 


2jt(f,n,  +  f2n2) 


+ 


Un„n3Sin| 


2jr  (f|n 


I  +  ^2*^2) 


(39) 


8.  Error  Analysis  Results 


This  section  contains  the  numerical  evaluation  of  the  cyclic  convolution  error  for  the  two  most 
prevalent  kernel  functions  in  geodesy.  These  kernels  are 


h(x„X2) 


(40) 


used  in  the  computation  of  the  geoid  undulation  from  gravity  anomalies  (planar  approximation  to 
Stokes’  function);  and 


h(X|,X2)=,  .  ..3,2 

(xj  +  x^) 


(41) 


used  in  the  computation  of  the  vertical  gradient  of  harmonic  functions  (upward/downward 
continuation  in  Molodensky’s  geodetic  boundary  value  problem);  it  also  appears  in  the  terrain 
correction  and  other  convolutions  (again,  always  in  planar  approximation).  To  simplify  the 
graphic  depictions,  the  computation  points  are  restricted  to  k|  =  k2  =  k;  the  frequencies  are 
restricted  to  f|  =  f2  =  f  ;  and,  N  |  =  N2  =  N  . 

Figures  3,  4,  and  5  correspond  to  the  geoid  undulation  kernel  (40);  while.  Figures  6  and  7 
correspond  to  the  vertical  gradient  kernel  (41).  Eaeh  figure  shows  the  magnitude  of  the  error 
transfer  function,  |  U).  |^(f,f)  | ,  depending  on  the  computation  point  (kAx,kAx)  and  the  frequency 
pair  (fj,fj) ,  where  fj  =  j/(2NAx) ,  and  N  =  96.  For  other  values  of  N  one  obtains  similar  plots, 
and  the  following  discussion  is  meant  to  be  more  qualitative  than  quantitative. 

Consider  Figure  3,  showing  the  transfer  function  for  the  cyclic  convolution  error,  (e])|^|^  • 
Each  vertical  (frequency)  profile  of  the  function  represents  the  transfer  function  of  the  error  for  a 
particular  computation  point  of  the  convolution.  Clearly,  they  are  like  low-pass  filters,  affecting 
primarily  the  lower  frequencies  of  the  signal.  Also,  as  the  computation  point  nears  the  edge  of  the 
signal  area,  the  cyclic  convolution  error  becomes  larger,  as  expected.  Nevertheless,  if  it  is  not  too 
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close  to  the  edge,  the  error  may  be  acceptable.  This  is  especially  true  for  the  kernel  (41)  (see 
Figure  6). 


0  4  8  12  16  20  24  28  32  36  40  44  47 

X  /  Ax 

Figure  3.  Transfer  function  of  cyclic  convolution  error  8 ,  for  kernel  (40). 

There  is  also  the  interesting  possibility  with  the  kernel  (40)  (Figure  3)  that  if  the  signal  has  no 
low-frequency  content  then  that  same  acceptable  error  occurs  for  computation  points  even  closer  to 
the  edge.  In  the  example  of  Figure  3,  if  x  =  12  Ax  represents  the  computation  point  at  which  the 
error  transfer  function  is  negligibly  small  for  all  frequencies  of  the  signal,  then  x  =  39  Ax  is  the 
computation  point  where  the  error  is  negligible  if  the  signal  has  no  frequency  content  below 
f=  5/(192  Ax). 

Figure  4  shows  the  transfer  function  for  the  cyclic  convolution  error  (£2)^15,  caused  by 
improperly  padding  the  kernel  function  (40)  with  zeros.  Such  a  procedure  has  the  potential  to 
reduce  computation  times,  compared  to  the  correct  padding  method;  but  as  seen  in  the  figure,  it 
should  only  be  used  if  the  signal  is  high-pass  filtered. 
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X  /  Ax 

Figure  4.  Transfer  function  of  cyclic  convolution  error  £2  for  kernel  (40). 

Figure  5  verifies  that  the  edge  effect  error  and  the  cyclic  convolution  error,  e,  ,  have  similar 
spectral  characteristics,  as  evidenced  by  their  transfer  functions.  Therefore,  in  avoiding  the  edge 
effect  by  restricting  the  computation  point  of  the  convolution  to  an  interior  subdomain  of  the  signal 
area,  one  also  avoids  the  cyclic  convolution  error,  even  if  no  zero-padding  is  performed. 

Figures  6  and  7  show  the  analogous  spectral  transfer  properties,  respectively,  of  the  circular 
convolution  and  edge  effect  errors  for  the  kernel  (41).  Because  this  kernel  is  much  narrower  near 
the  origin,  the  errors  are  strongest  only  for  computation  points  very  close  to  the  edge  of  the  signal 
area.  Note  the  break  in  the  abscissa  scale  in  these  figures  and  the  similarity  in  the  transfer 
functions. 
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Figure  5.  Transfer  function  of  edge  effect  error  for  kernel  (40). 


Figure  6.  Transfer  function  of  cyclic  convolution  error  8}  for  kernel  (41). 
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Figure  7.  Transfer  function  of  edge  effect  error  for  kernel  (41). 


9.  Numerical  Example 

The  circular  convolution  error,  E| ,  given  by  (32),  and  the  edge  effect,  Eejjgc,  given  by  (34), 
were  calculated  from  a  regular  grid  of  equal-area  mean  gravity  anomalies  convolved  with  the  kernel 
(41 ),  This  convolution  yields  the  vertical  gravity  anomaly  gradient,  given  in  planar  approximation 
by 


oo 


Ag(x',y)-Ag(x,y)  x^dy- 

((x-xf +  (y-yf ) 


(42) 


The  total  data  array  is  a  grid  of  2  N  j  =  680  by  2  N2  -  540  values  in  the  west-midwest  part  of 
the  United  Slates.  The  grid  mesh  size  is  4  km  by  4  km.  The  circular  convolution  errors  and  edge 
effects  were  computed  for  the  inner  N  j  x  N2  area  and  are  shown  in  Figures  8  and  9,  respectively. 
As  predicted  by  the  error  transfer  functions  (Figures  6  and  7),  the  errors  are  largest  near  the  edges 
of  this  computation  area.  The  most  inner  N  j/2  x  N2/2  area  shown  in  Figures  8  and  9  by  dotted 
lines  indicates  the  area  that  might  be  considered  free  of  edge  effects  (although  it  could  be  larger). 


-  17- 


In  Table  3,  the  root-mean-square  (rms)  errors  are  listed  for  each  type  of  error  and  for  each  region, 
inside  and  outside  the  dotted  lines.  It  is  clear  from  the  figures  and  the  table  that  the  circular 
convolution  error  and  the  edge  effect  are  commensurate.  Avoiding  the  latter,  which  must  be  done 
in  any  case,  also  avoids  the  former,  and  no  zero-padding  of  the  data  is  needed. 


Figure  8.  Circular  convolution  errors  (absolute  values)  for  the  vertical  gradient  of  the 
gravity  anomaly.  The  contours  decrease  in  value  from  the  edges  with  interval  equal  to  0.2 
Eotvos. 
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Figure  9.  Edge  effect  (absolute  values)  for  the  vertical  gradient  of  the  gravity  anomaly. 
The  contours  decrease  in  value  from  the  edges  with  interval  equal  to  0.2  Eotvos. 


Table  3:  RMS  values  of  the  circular  convolution  error,  edge  effect,  and  their  sum  for  the  example 
of  the  vertical  gravity  anomaly  gradient  described  in  the  text.  Units  are  Eotvos.  Areas  are 
delineated  by  the  dotted  lines  in  Figures  8  and  9. 


Edge  Effect 

Circ.  Conv.  Err. 

Total 

RMS  of  Inner  Area 

0.023 

0.027 

0.034 

RMS  of  Outer  Area 

4.478 

4.434 

6.116 

Abs.  Max.  Err. 

135.2 

129.4 

167.3 

10.  Derivatives 

Derivatives  of  signals,  and  of  convolutions,  play  an  important  role  in  geodesy.  For  example, 
consider  the  deflection  of  the  vertical  according  to  the  Vening-Meinesz  formula,  that  is  the 


derivative  of  Stokes’  formula,  that,  in  turn,  is  a  convolution  of  gravity  anomalies  and  Stokes’ 
function.  As  another  (rather  arcane)  example,  the  terrain  correction  is  the  vertical  gravitational 
attraction  at  a  point  P  due  to  the  mass  attraction  of  the  residual  topography,  Ab  =  b  p  -  bg  ,  with 
respect  to  point  P.  In  planar  approximation,  it  is  given  by  a  series  of  iterated  convolutions 
(Sideris,  1990).  Klose  and  Ilk  (1993)  give  the  correct  formulation  of  the  attraction  in  terms  of 
Fourier  transforms; 


A(P)  =  kp  i  A,(P) 
1 


(43) 


where  k  is  Newton’s  gravitational  constant,  p  is  the  (constant)  density  of  topographic  mass,  and 


A,(P) 


(^2r-l 
(2r)!  n  =  0 


(44) 


where  F  (  • )  denotes  continuous  Fourier  transform,  and  f  is  frequency  given  by 

f  =  V^fl  +  f5  (45) 

The  factors  f^'^'  '  have  a  conventional  (though  not  required)  interpretation  (Sideris,  1990)  as  being 
due  to  vertical  derivative  operators  applied  to  the  residual  topography,  formally  extended 
harmonically  into  exterior  space  like  a  potential. 

Finally,  having  computed  the  terrain  effect  according  to  (43)  (or  some  approximate 
convolution),  it  is  possible  using  standard  properties  of  the  Fourier  transform  to  calculate 
derivatives  of  the  attraction,  thus  yielding  terrain  effects  for  the  gravitational  gradients.  It  is  first 
noted  that  from  the  relationship  between  A(P)  and  the  corresponding  generating  potential,  T,  that 
is,  A(P)  =  -f')T(P)/3z  ,  we  obtain 

F(T)  =  ^F(A)  (46) 


Then,  the  second-order  gradients  of  the  potential  can  be  expressed  formally  as  (see  also  Tziavos  et 
al.,  19SS) 


r  -  _p-i 


(i27ifn)(i27tf,,) 


2jr  f 


F(A)  ;  n  =  1,2  ;  m  =  1,2 


=  '((i27ifjF(A));  n=l,2 


dx  „  dz 


(47) 

(48) 
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(49) 


In  each  of  the  cases  listed,  (44),  (47),  (48),  (49),  and  others,  the  desired  quantity  is  the  inverse 
transform  of  the  product  of  a  two  functions  of  frequency,  one  representing  a  signal  and  the  other 
representing  derivative  operators.  Thus,  according  to  the  convolution  theorem,  the  desired 
quantity  is  a  convolution  -  but  only  formally  in  these  cases  since  the  derivative  operators  must  then 
be  viewed  as  derivatives  of  the  Dirac  delta  function.  When  approximating  these  continuous 
“convolutions”  with  cyclic,  discrete  convolutions  (amenable  to  FFT  processing),  where, 
furthermore,  F(A)  may  be  already  approximated  by  a  cyclic,  discrete  convolution,  it  is  important 
to  define  properly  the  discrete  versions  of  the  derivative  operator  “transforms”,  such  as  (i  27t  f„)  in 
equation  (48). 

Since  the  signal  being  differentiated  is  assumed  to  be  real,  the  result  of  the  cyclic,  discrete 
convolution  with  the  derivative  operator  should  also  be  real.  Therefore,  the  DFT  of  the  signal  and 
the  DFT  of  this  convolution  both  satisfy  the  conjugate  symmetry  property  (cf  (5)),  which  means 
the  DFT  of  the  derivative  operator  should  also  satisfy  the  conjugate  symmetry  property.  For  two 
dimensional  periodic  discrete  sequences  defined  over  the  fundamental  area, 

-N,/2<(>,<N,/2-1  ,  -N2/2<(2<N2/2-1  (50) 

the  conjugate  symmetry  property  is 


(51) 


The  necessary  condition  (51)  and  periodicity  of  the  discrete  spectrum  imply  that  four  spectral 
values  must  be  real: 


(^n,.nJo,o  “  (^n..nJo_o  ’  (^N„Nj„__Ny2 " 


N:/2 


2,0’  (^N„nJ_n 


i/2 -N 2/2  i  ^  ^  1  2  /-N  ,/2,-N  2/2 


In  addition,  there  is  conjugate  symmetry  in  l’|  for  C2  =  0  and  IS  “  ”^2/2  • 

(^N,.N  0  "  (^N,.N  "  (^N,.N2S_N3/2 

and  in  (2  for  ('1  =  0  and  (|  =  -N 1/2  : 


(52) 


(53) 
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(^Ni.nJop^  (^N„nJ(,__(i^  >  (^N,,Nj_N,/2,f2  (^N,.nJ_n |/2  _ 


(54) 


Thus,  given  the  discretization  of  the  derivative  operator  transform  for  the  frequencies  (52),  one 
must  ensure  that  symmetries  (51),  (52),  (53),  and  (54)  are  fulfilled.  This  is  not  automatic,  for 
example,  in  the  case  of  (i  2:r  f„) ,  and  one  must  override  the  continuous  definition  with  the  discrete 
symmetries.  Then,  for  FFT  processing,  the  resulting  transform  may  be  translated  by  periodicity  to 
the  domain  0<(|<N|-1,  0<C2^N2-1  using  a  procedure  analogous  to  (28). 
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Fig.  10.  Discrete  kernel  of  the  derivative  operator  3V(3x|  8x2)  (N 1  =  N2  =  32  ). 

As  one  might  expect,  the  cyclic  convolution  errors  are  significant  only  for  computation  points 
close  to  the  edge  of  the  computation  area.  Furthermore,  because  of  the  ripple  effect,  the  error 
dominates  in  the  high  frequency  domain.  Again,  to  avoid  the  cyclic  convolution  error  one  should 
extend  the  signal  with  zeros,  as  in  (25)  or  (27).  On  the  other  hand,  the  edge  effect  remains,  and 
will  be  of  the  same  order  of  magnitude  as  the  cyclic  convolution  error  if  the  signal  is  not  extended. 
If  the  derivative  operator  transform  multiplies  a  convolution  transform,  as  in  (47)  -  (49),  and  it  is  a 
cyclic  convolution  operating  on  extended  (padded)  arrays,  then  the  derivative  transform  should  be 
extended,  as  well,  and  should  be  padded  with  its  own  natural  values  according  to  the  symmetry 
properties  (51)  -  (54). 
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X  /  Ax 

Fig.  11.  Transfer  function  of  cyclic  convolution  error  Ej  for  the  kernel  of  Figure  10. 


11.  Summary 

The  methods  of  approximating  a  convolution  by  a  discrete,  cyclic  convolution  (for  FFT 
implementation)  in  the  past  have  been  studied  empirically  as  applied  to  geodetic  problems.  Various 
techniques  were  tested  on  actual  data  by  numerous  investigators  to  determine,  in  some  probabilistic 
sense,  the  accuracy  of  the  methods. 

The  results  of  the  present  discussion  are  three-fold.  Expressions  are  derived  for  the  discrete 
circular  convolution  error,  for  the  error  committed  by  improper  zero-padding  of  the  kernel 
sequence,  and  for  the  edge  effect  (discrete  approximation).  Second,  through  the  use  of  error 
transfer  functions  and  a  numerical  example,  it  is  shown  that  the  cyclic  convolution  error  committed 
by  not  extending  the  arrays  is  not  greater  than  the  error  associated  with  the  edge  effect.  Therefore, 
by  limiting  the  computational  area  to  reduce  the  edge  effect  (that  does  not  disappear  with 
extensions),  one  automatically  reduces  by  similar  amount  the  cyclic  convolution  error.  This  is 
important  when  convolving  large  arrays,  and  extensions  (quadrupling  memory  requirements)  are 
not  feasible.  A  by-product  of  the  analysis  is  that  extending  the  kernel  with  zeros  makes  no  sense 
from  an  accuracy  viewpoint,  although  the  error  is  less  than  with  no  extensions.  And  third,  the 
analysis  is  extended  to  properly  defined  discrete  derivative  operator  transforms  with  a 
demonstration  through  error  transfer  functions  that  the  resulting  cyclic  convolution  error  is 
confined  to  the  high  frequencies  and  to  computation  points  very  close  to  the  edge,  as  expected. 
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